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1. Introduction 


The numerical solution of the initial value problem for u system of 
ordinary differential equations 

y* (t) - f <y(t),t) (la) 


y(t o ) 


o 


(lb) 


can be considered locally as the computation of a sequence of approximations 

Y ■ (y , .,,...,y . ,y } with each y beinq a numerical approximation to 

n n-ktl 'n-l 'n n 

n 

y(t ) for t ■ t + F. h.. While many sophisticated packages f 1,2# 3. 4 , 
n n o , , i 

i=l 

exist which change the stepsize h. in order to achieve stability and specified 

accuracy, the stepsize usually remains constant for a number of st^ps, and tin 

the change of stepsize is usually accompanied by an intcr|<olatory change in 

the solution sequence Y , although not always [5 j. For these reasons the 

n 

local analysis, limited to k steps, given in the sequel will assume constant 
stepsize h. The actual numerical solution element y^ is usually computed 
from elements of the sequence Y^ ^ and f(y,t) using oni of several formulas, 

Jf -f 1 

each of which lias a local discretization error term <t>^Ch where <p depends 

on f aid L , C is a constant dependent, on the formula, and r is called the 
n 

order of the method. This work is concerned wi h how errors are propagated 
when numerical methods are applied to a differential function f(y,t) which 
is nonlinear in y, witn emphasis on stable propagation of the errors. 

The most typical formulas in use are Runge- Kutta multistage formulas 
[f>] and Adams type multistep methods. An s-st.ige Runge-Kutta formula is 

(2a) 


K * hf(y ,t ) 
o n n 


q-1 


< 1-1 


K - hf(y + F B .K., t + ' h£ . ), q * 1,2 s-i 

q n . „ gj j n . . qj n 

j=0 3=0 

s-1 

y .. * y + E Y K . 
n+1 'n q q 

q-0 


(2b) 


(2c) 


families of explicit formulas exist for orders one through four, with 
thi order equal to the number of stages. Fifth order formulas require at 


leas* - six stages, ami more accurate implicit formulas formed b/ replacing 
(2b) by 

‘1 * M * 

K » hf(y + E 8 . K., t + T. hP . ), and requiring solution of 

« n J .. 1 n , 0 

system;; of nonlinear equations, also exist [7]. The Adams-Bashf<<rth and Adam:.- 


Moulton formulas are of the form 

k 

y - y + hJi b. f (y . ,t ) , 
n n-1 . , i n-i n-r 

i~l 


y = y , + hJ, b. f (y . ,t ) 

n n-1 i n-i n-x 

i=0 


( 3 ) 


(4) 


respectively. A k-step Adams Bashforth formula is of order k, and a k-step 
Adams-Moulton of order k+1. Since the latter may require iterative solution 
of a system of nonlinear equations, the initial solution guess may bo pre- 
dicted by an Adams-Bashforth formulas of the same or one lower order, and 
the difference may be used to approximate the discretization error [ 6 J . 

Many methods iterate (4) a set number of times, often one, rather than to 
convergence. 

Th'! concept of multistep methods is useful since order of accuracy as 

high as k+1 [si can be achieved by a stable formula with only a few 

evaluations of f(y,t), whereas a correspond i ugly accurate kunqe-Kutta 

formula requires at least one evaluation of f(y,t) for each additional 

order of accuracy. Newer methods have been developed to handle special 

circumstances. The backward differentiation implicit formulas 

k 

y = E a. y + h b f(y ,t ) (5) 

n , l n-1 o n n 

x=l 

are useful in solving stiff systems of equations [1]. For systems witli 
known complex time constants \, such as are often encountered Ln simulation 


of physical py stems such as aircraft, formulas dependent on h\ toi constant 

h arc known. Such formulas have the general form 
k k 

»n ’ f Vl ♦ hr (f ’> 

i"0 

The coef f iciiMit s •• | *h | are chosen to lit exactly tin? case y(t) 

I'^e* coayt + sinyt where A-H'ly and It Is constant. The order of these 

methods is generally two lower than would l>e expected since two degrees of 
freedom are lost in forcing the solution to follow $(t). 


The above indicates some of the wide variety of problems which can he 
handled by multistep formulas. However, a better understanding of the 
stability problems inherent in using them would be helpful. 


2. Stability of Solution Sequences 

The standard linear stability analysis is formula specific and 

describes, the behavior of a numerical formula applied to the complex test 

equation y' - Ay , y(t ) f 0. Let B ■» {e(t ),..., e(t , ) , e (t ) } 

o n n-Kri n-1 n 

bo the difference sequence tor y(t )-S(t ) where each solves the test 

n n 

equation for a different initial value y . z . Then e(t ) ® (y )e' 1 

Do n 0 0 

is nonincreasing in norm for Re (A)<0. Such a condition is called stability. 

It is desirable for the numerical solution y to bo stable if the true 

it 

filiation is stable, so for a given stepsize It, one finds all complex A such 

that any numerical sequence K ■ {e , , .,<> } has the property 

it n-k+l n-1 n 

| cv | for .ill i, where e ^ *= yv-z., the tlifforence between 1 1 to numeri- 
cal solution sequences with initial values V 0 ».' ( ,- 

For ruler's formula, (l4h\)y , hence the numerical solution it. 

stable for |l+hA|<l. lor multistep fornutlas, linear stability is character- 
ized by the genoratinq polynomials 


<7a) 


k 

P(C) ■ X* 

f k-i 

a. f. 

i-0 

k 

0(C) - T. 

i rt 

. r k_i 
b i L 


where p and 0 have? no common diviners. The? region of linear stability 
is all hX such that p(C) * hA o(C) lias all roots inside the unit circle, 
or on the unit circle and simple [6|. For Filler's formula, p(C) ♦ hX o(C) 
is (-f,+ l) *■ hX. Since this analysis is formula specific, tc investigate 
the formula's effect on an actual f (y,t) one considers all the eigenvalues 
X. of the Jacobian matrix <*f (y, t)/3y j if all hX ( are inside the stability 
region for all y^, t of interest, then the numerical solution will be 
stable. Insuring tliat such a condition holds is usually not desirable and 
often is impossible. 

To develop a stability analysis for nonlinear f(y,t), let f(y,t) have 
the property 

Re <y-z, f (y, t) - f (*,t)> < n| | y-*| \ ? (H) 

for all t, y, and /. of interest. Here <u,v> » U*Qv for some positive 

•j 

definite Hurmitian matrix Q, and || u | | - <u,u>. Then for any two 

solutions 


y(t), z(t) ■ y(t)-e(t), e(t) satisfies 


- f (y(t) ,t) - f (z(t) ,t) 


and (0) implies that 


dt 


||e(t)||*' - 2 Re<e (t) , < 2u||e(t)|| 


and thus 


H«(tl|| £ e Ut ||e(t o )||, 
which is non-increasing for H<0. 

However, (H) is again a condition that cannot be easily verified, so 

a concept relating the true solution sequence Y* = { y ( t 

n n-k-t 1 


y ( t - n _ 1 ) » y lt n ^ to tl,e commuted sequence Y^ - (y^_^ f ^ , . . .y^ ^ ,y^ } is 


presented in the sequel. The following definitions and theorem are 


helpful. They occur in Dahlquist 10) and the theorem proof is presented 


in Brown 110). 


Definition: a linear k-step formula satisfies 


°-;. 0 u )V) * •> v'vj'Vi 11 


Definition: a one- leg k-step formula corresponding to (9) satisfies 


k k k 

0 = E a.y . + lisf(— E b.y . , — E b,t .) 

j-o J n “J s j-0 j n * j s j-o i n -i 


where s ** o(l). Without loss of generality, set 8=1, 


Theorem 1 [9] Let Y be a sequence which satisfies (10), and let 


Y = {y } be such that 
n n 


y = E b.y . » o(E) y , 
n jM0 j n-j r n 


where E denotes the back shifting operator. Then Y satisfies (9). 

n 


Conversely, if Y satisfies (9) then there exists a sequence Y such 
n n 


that y = o(E)y , and y satisfies (10). 
n n ' n 


This shows that Y given by the one-leg k-step formula will 
n 


have similar stability properties to its corresponding linear k-stop 


sequence Y. Dahlquist (9) has described a discrete Liapunov function 


V„ p i which, applied to a sequence Y , characterizes the stability 
Ci,c,n n 


i i 


of that sequence generated by a non-linear system (1). Let 

t t 

V G,t,h <V ' ^ °ij <V n-l+l* Vj.l’ whcr ” 0 18 8 

positive definite, l x t matrix. The structure of G assures that 

V , is positive for Y / { 0 } . 

G,£,h r n 

Definitions The (G,l h) - don.) in of attract ion of (the numerical 

solution to) the system (1) is all z such that AV p (Z ) « 

0 t*,t,h 0 

V G,£,h ( V " V G,£,h ( V- 0 ' Where Z 0 = .* 0 }. ^ 

{z( (2— £) h) , . . . ,z Qt Zj } for the numerical solution and ^ = (z((2-l)h), 
z(h)} for the exact solution. 

Definition: The (G, £, h)-stability region of (the numerical 

solution to) the non-linear system (1) is all z Q such that 

V p . (Z„) < inf (V„ » .(Z-) } where 3D is the boundary of the 
G,t,h 0 - G,£,h 0 

*o e,D 


(G,£,h)-domain of attraction. 

This has the following application. Rather than requiring that 
Re<y-z, f (t,y ) -f (t . z) > <_ 0, a connected subset of initial values y () . 
found such that y(h) will be in that subset if y Q was; this is the 
stability region. This insures that the difference y(h)-z(h) is 
bounded since both y (h) and z(h) are in the stability region if y Q 


and z„ were. If f(y,t) is autonomous, y(t ) will remain in the region 
0 n 

as n ►««’. For most well behaved functions f(y,t), the boundary of the 
region around a stable point can be approximated computationally. 

Once the analytic stability region is known, the numerical stability 
region can be calculated using the one-leg k-step method for the same 


The two regions can 


sequence {y (1-k) h) , . . . ,y (-h) ,y Q ) to get y^. 
then be compared. 

Analytically, it is possible to form a particular G, based on 
the coefficients of a one-leg k-step method, such that all numerical 
sequences based on f(t,y) that satisfy (8) will have a stable solution. 
Liniger and Odeh [ll] have shown how to pick G for second order two-step 
formulas, second order three-step formulas, and third order, three-step 
formulas. 

It is shown below that even or. arbitrary choice of the positive 
definite hermitian matrix G will generate some usable results, and 
theorem 2 demonstrates that using some G for a one-leg k-step solution 
will generate the same stability region for the related solution 

of the linear k-step formula for a modified G. The proof appears in 
Drown (10], 

Theorem 2 If V. „ , (Y ) * c for the rymmetric positive definite 
c,n n 

A 

matrix G, then there exists a symmetric, positive def in ite matrix G, 

dependent only on G and o(x), such that V* ^ j/'n ** c ' w * iere 

Y » o(E) Y are the elements of Y and Y . 
n n n n 

3. Existence of Stability Regions 

Sufficient conditions can bo developed for the existnece of the 
(G,l,h) -stability regions based on known techniques such as Liapunov's 
direct method l 121 and property (8) . An important concept in the 

A 

development is the equilibrium y (t) of the differential function f(y,t). 

A 

While some references define it for an initial value y(t Q ) ■= 0 in the 
space of the dependent variables, this is accomplished by an unnecessary 


change of variables that could be confusing. The important point 

A A 

is that y(t) satisfies f(y(t Q ), t Q ) ■ 0, so that, if f is autonomous, 

then y(t) is a constant, and otherwise, the Taylor scries about t () 

* 2 

is given by y(t C) ) ♦ (t-t Q ) f ' (*)/2 and is thus slowly varying and 
nearly constant for t near t^. 

AAA A A 

Definition: The solution y(t) of y' ■ f(y,t) for y(t Q ) ■ y Q , 

A 

such that f (y Q ,t 0 ) ■ 0, is called the equilibrium of f(y,t) . 

Definition: The equilibrium of f(y,t) is said to be asymptotically 

stable if there exists a ctt^,®) such that for every e>0 there 

A A 

exists = 6j(e,T 1 )>0 such '.hat if l^o~ y o^ <4 l' then I I V "V <t) | I <c 
in (t^, °°) , and there exists a x.,e (t^®) and 6 ,(x,) such that if 

A A 

| J y 0 — Yo'l I <fi 2 2 ^ then lim I |y (t)-y ( fc ) I I ■ °* Where y(t) satisfies y * ( t ) 

t-«» 

- f(y,t), y(t Q ) » y Q . 

Theorem 3 (Liapunov [12J) The equilibrium of f(y,t) is asymptotically 
stable if there exists a function v(y,t ) which is positive definite 

A 

in some region D about y(t Q ) and lim v(y,t) = 0 uniformly in t as 

A 

| | y-y | | -K) there, and whose total derivative is negative definite on D. 

With this background, it can be shown that a well behaved f(y,t) 
which has a not necessarily unique Liapunov function v(y,t) with region 
D implies the existence of a ( I, t,h ) -domain of attraction and stability 
region of both the exact solution and of a one-let k-step solution 
based on a stable multi-step formula. Similar proofs can be developed 
fro positive definite matrices G ^ I. 

A 

Theorem 4: If the equilibrium y(t) of f(y,t) has a Liapunov function 

of the form v(y,t) = y*Qy for a positive definite matrix Q, y* the 






u T 

• r 


transpose of y, and f(y,t) is continuously differentiable on a convex 
domain D whose boundary 3D is defined by v'iy^jt^) 0, and il f ( y , t > 
has the property on D that ty j -y ,) If iVj ,t ) -f (y , , t ) ) mi | |y J ■ y^ 1 1 ", 
where ii<0 and x*yjx ■ ||x}|‘, thrn for any ixsint y^. • y(t^) in the 
interior of P, an (I ,( , h)-stabili ty region can be constructed using 
rays for y^, provided h <_ h^if.f). 

P» oof Note that the solution from to t ( ♦ h is spiraling in from 
30 toward the equilibrium, which is inside a circle | |y(t)-y ||, • 

•% U • “ 

— M where M * max f'U), since the hypotheses and (8) gives 

«D 

llytt^) ” | li exp(ufh) I |y(t 0 )-y(t Q ) 1 1. By dofiniticn, AVj ( h (Y^) 

■ I lyt^ll^- I |y(t Q ) | |‘, which occurs when 0 - viy^., t^)-y(y 0> t 0 ) 

■ (y^*y,)v'iz) by the mean value theorem. Since v* (z) 0 on the 

boundary 3D, along any ray from y ( . one can fin*< yit^) generated by a 

y Q or. 3p, provided hf < is small enm ih that the trajectory from 

<v <>, 

y 0 to y^. is entirely in D. 

The boundary of (I,f’,h) -domain of attraction D* can be constructed 

of all such rays. The (I,l,h)-stability region P‘* has a boundary of 

all points such that V , » v* » min V, . , , which can also be 

I,t,h 3D , l,f,h 


constructed. 


3 d* f(* ) 


Theorem 5t It the hypotheses of theorem 4 are met and 1 — A-‘ — 

3y dt*^ 

is continuous in D, then for any y ( . in the interior of P, an (I, (, In- 
stability region can be constructed for a p-th ordei one-leg k.-stej 
method, for h <_ if , f ) . 


rrw - f - AV - ||y f ||* - | |y(t f -Ch) | I ’’ where 


1 L 




4 l A. 




k k k 

y / m l a 4V +hf( l Vlt/.J* l b 

L j-l 3 43 j-o 3 c 3 j-o - 4 3 

and y(t^) ■ y^, ♦ h I>+ ^f^^(2)» the truncation error formula. 

If the truncation error varies continuously as y(t 0 > varies along 3D, 
then two solutions and z^. generated by y 0# z 0 on 3D have the 
property that y^-* z ^ uniformly as y Q -+ z Q , and a smooth curve 3D' exists 
on the boundary of the (I . i, h, ) -domain of attraction. Similar arguments 
show the existence of the (I ,£,h) -stability region. 


•1. Example 


The system [12J 
x ' ■ -x 
y' ■ 5x 


2y + x“ sin(t) 
y + y/(t+4) 


( 12 ) 


has a Liapunov function for Q ■ q.^, ^ B 37/44, q,., ■ 4/11, 

q ^2 ■ q, j ■ 3/44. Figure 1 shows D, as well as D', the (I, 2,0.1)- 
domain and D" , the (1 ,2, 0. 1 ) -stability region. The solution was 

CD 

generated by assuming a power series expansi- n x(t) ” £ a^t 3 ,y(t) ■ 


j“0 


} b.t J for which a - x , b Q ■ y~, and a. and b. can be generated 
j-0 3 u u u u 1 * 

recursively. The series converges for | t 1 1. Figure 2 shows D' and D" 
for the numerical method based on the Adams-B.ishforth 2-stop formula 


y n+l “ y n + h<1.5*f(y n ,t n ) “ * 5#f <Vl' t n -l ) ) * 


(13) 


REFERENCES 


1. Goar. C. W. , Algorithm 407s DIKSl'D for Solution of Ordinary Differential 
Equations, CACM 14, pp. 185-190 (1971). 

2. Hindmarsh, A. C., GEARs Ordinary Differential Equation System Solver. 
UC1D-300001, Rev. 2, Lawrence Livermore Lab., Livermore, California (1972). 

3. Krogh, F. T. , A Variable Step Variable Order Multistep Method for the 
Numerical Solution of Ordinary Differential Equations. In Inform ation 
Processing 68, Vol. I, pp. 194-199, A.J. H. Morrell ed., North Holland, 
Amsterdam (1969). 

4. Shampine, L. F. and Gordon, M. K., Computer Solution oi Ordinary Dif^r- 

ential Equations: Initial Value Problems, Freeman, San Francisco (1976). 

5. Tu, K. W. , Stability and Convergence of General Multistep and Multivalue 
Methods with Variable Stepsize, UIUCDS-R-72-526, U. of Illinois, Urbana, 
ILL. (1972). 

6. Gear, C. W. , Numerical Initial Value Problems in Ordinary Differential 
Equations, Prentice-Hall, Englewood Cliffs, N J (1971). 

7. Cash, J.R. , Semi-Implicit Kunge-Kutta Procedures with Error Estimates for 
the Numerical Integration of Ordinary Differential Equations, Journal ACM 
23, pp. 455-560 (1976). 

8. Dahlquist, G., Numerical Integration of Ordinary Differential Equations, 
Journal ACM 23, pp. 455-5bO (1976). 

9. Dahlquist, G. , On Stability and Error Analysis for Stiff Non-Linear Pro- 
blems, Report NA-7508 , Dept, ot Information Processing, Royal Institute 
of Technology, Stockholm (1975). 


’0. Blown, R.L., Stability of Sequences Generated by Nonlinear 
Differential Systems, to appear* Math. Comp. 

11. Linigor, W. and Odeh, r., On Liapunov Stability of Stiff 
Non-linear Multi-step Difference Equations* AFOSR-TR-76-1023, 
IBM Thomas J. Watson Research Center (1976). 

12. Lehnigk* S. H., Stability Theorems for Linear Motions, Prenti 
Hall* Englewood Cliffs, N.J. (1^06) 





